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ABSTRACT 


In this paper we consider acoustic shock waves in a variable area duct 
which contains near sonic flows. The problem we treat here is modelled 
after an aeroengine inlet. It is known experimentally that area variation 
of a duct and high Mach number mean flow can reduce acoustical energy 
yielding substantial noise reduction. One possible reason for this is acoustic 
shocks. We describe the use of an explicit numerical method which is very 
accurate and also captures shocks reasonably well. Comparisons of the 
results are made with an existing asymptotic theory for Mach numbers 
close to unity. When shock occurs reduction of sound pressure levels are 
shown by examples. 
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1. Introduction 


Recently a numerical solution for the propagation of sound in a vari- 
able area duct which contains a high Mach number subsonic flow was 
studied by the authors (reference 4). The nature of the wave propagation 
was nonlinear. This note is a continuation of this work and presents some 
results concerning computation of shocked waves in this situation. The 
model problem studied here and in reference 4 serves the purpose for the 
study of an aero-engine inlet. Briefly there is a flow in a variable area duct 
and acoustic waves propagate upstream of this flow. The acoustic waves 
are generated by an incident plane wave on the left of the duct and leave the 
right end without reflecting. It was first observed experimentally (reference 
3) that with a given flow and a proper choice of area variation it is possible 
to attenuate the sound intensity as much as 20db. Theoretical reasons for 
this mechanism of reduction of noise level are still under investigation. 
They are complicated by the fluid equations which are the Navier-Stokes 
equations. However, in the inviscid limit it is believed (reference 7, 8) 
that acoustic shocks which result in energy loss is one such possibility. 
Though sound level attenuation was not shown in ref. T, they did show 
energy loss when shocks occur. The work in reference T is based on an 
asymptotic theory and valid for Mach numbers close to unity. In a later 
experimental study (reference 5) it was shown that even for Mach numbers 
far less than unity (about 0.7) one still obtains substantial sound reduction. 
Thus a numerical study was undertaken to calculate solutions for all Mach 
numbers. 

There has been only a little published in the literature for this flow 
configuration. The classical Fubini solution uses an asymptotic expansion 


to obtain an approximate solution to the one dimensional gas dynamic 
equation for a uniform duct. Polyakova (reference 10) made extensions for 
problems with flow and Blackstock (reference 1) for the case of shocks. For 
variable area ducts Myers and Calleghari (reference 8) used the method of 
matched asymptotic expansions. This reference was the only source for any 
constructive shock solution in the literature. Recently parallel to our work 
Walkington and Eversman (reference 12) carried out computations of this 
situation. Our study yields similar results, but we believe, our approach is 
simpler for the reasons indicated in the next section. Moreover the scheme 
used here is more accurate than the one in reference 12. The authors 
would like to point out that there has been other work namely Nayfeh. et 
al (reference 9) to compute the nonlinear, but unshocked solutions for this 
type of configuration. 

In this work we describe the model briefly and indicate the assembly of 
the numerical method. We emphasize the approach we took for obtaining 
boundary conditions. We present numerical results obtained for shock 
cases and compare with the asymptotic results of reference 7. We also 
present noise level distribution over the duct and demonstrate the noise 
reduction in the situation of a shock. 

2. Equations of motion 

The total flow field is governed by the inviscid, compressible Euler 
equations. They consist of equations for continuity, momentum and 
energy. The energy equation can be replaced by isentropic relations 
provided one is seeking only weak shocks. This is exactly the case in 
acoustics where only weak shocks are the central goal. Strong shocks 
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cause disturbances in the main stream and the meaning of acoustics will 
not be valid. This philosophy was adapted in reference 4 and it yielded a 
system of two equations for acoustic density and velocity rather than three 
equations thus reducing computational costs. 

The situation which is of interest here assumes a quasi one dimen- 
sional flow. The flow configuration is depicted in figure 1. In this a steady 
flow is moving from the left to right and the acoustic waves are propagat- 
ing up-stream of the flow from a harmonically varying source (plane wave). 
Then the total field is governed by the following equations, where A(x) is 
the area variation of the duct: 



ip u ) + 


p u dA 

A dx 


du dfu 2 7 

m + sAj + 7=1^" J 


= 0 


and the pressure is determined from 


( 2 . 1 ) 

(2.2) 


_ P o_, 
P = ~^P 
Pi 


(2.3) 


Here ~p,U and p are the total density, velocity and pressure, respectively, 
and the quantities with subscript zero denote ambient values. We divide 
these flow quantities into ‘mean’ and ‘acoustic’ parts. That is, if mean flow 
quantities are assigned a subscript s then: 

U — Ug + U 

P = Ps + P (2.4) 

V = Ps + P 

The mean flow is assumed to be steady and they satisfy the following 

steady state flow equations: 

d . . , p e u s dA 

5^rf.) + — ^ = o 
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(2.5) 


d_ 

dx 


( 2 . 6 ) 


U‘ t 


.[ _JL -j- -J—Elpi-l 


0 


2 7 - 1 PT 6 

V* = (2.7) 

Then equations (2.1) - (2.3) yield the following nondimensional acoustic 


equations. 


Pt + ( Usp + p s u + Up) x + ( u s p + p s u -f- up) 


l_dA 
A dx 


0 ( 2 . 8 ) 


u t + 


u 


< u + Y + 


P/ Pe + 


7-2 


p = ci 1 + 


(■ 


7~ 


(p//>s) 2 )] =0 (2.9) 

-p/p c )p (2.10) 


Here c 6 is the local sound speed in the flow and is given by 



( 2 . 11 ) 


The details of derivation of these equations are available in reference 
4. Note that in these equations density, velocity are scaled by their ambient 
values pQ, Co , the pressure is scaled by PoC'o and fc he area by the throat 
area At of the duct. Moreover the distance and time are scaled by and 
by W , where 0J is the frequency of the source. This frequency corresponds 
to the engine noise source. In this nondimensionalization process the 
meanflow speed becomes the Mach number distribution and the solution 
space becomes the interval [0,L], where L is the duct length multiplied by 

ft- 

These equations are to be solved subject to the following boundary 
conditions. At x = 0, 


«(0, «) = /(<) (OO). (2.12) 
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B{u,p)(L,t) = 0 


At x = L, 


( 2 . 13 ) 


where f(t) is the source variation and (2.13) dictates a nonreflective 
(impedance type) condition. We shall derive this condition in section 4. 
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3. Numerical scheme 

We shall briefly indicate the numerical method used in this work. 
This method is an extended version of MacCormack’s method with fourth 
order accuracy in space and second order accuracy in time (see reference 
2). However, it should be viewed as a fourth order accurate scheme, 
since we are interested in steady time averaged quantities. Let the spatial 
discretization of the axis of the duct be given by Xj = 

(j — 1 )L/J (j = 1 Let us define forward and backward flux 
difference operators by 

PfU) = Vj - 8/,±, + f j± 2 (3.1) 

Here the plus ‘sign’ denotes forward and ‘minus’ sign denotes backward 
operations, respectively. Then for a single equation of the form 

Ut -f* fx — h (3-2) 


the scheme works as follows: 


xfp — Uj 


f n ) _i_ ^ —h n 

6Ax Pj (/ 2 


«? +1/2 = \ k + + f^ + (/ ( 1 ) ) + t*?’ 


(3.3) 


This has a backward predictor and a forward corrector. In the next At/2 
time step it is changed into a forward predictor and a backward corrector 
as follows: 


uW = «y+>/ 2 + 

■»+• = I [af+'Z 2 + „(*) — ^Zip - ( /(*)) — 

2 2 p + 2 6Az 1 V ) ^ 2 



(3.4) 
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In these formulas subscript n,n - j— 1/2 and n \ denote quantities 
evaluated at times nAt,(n + 1/2)A t and (n + 1)A t and the superscript 
(1) denotes the predicted values. It is pointed out in reference 9 that 
alternating formulas (3.4) and (3.5) at each time step is necessary to achieve 
fourth order accuracy for nonlinear problems. 

We note in (3.3) and (3.4) the flux difference operator P^~ is not 
defined for j = J-l and j= J and P~ is not defined for j=l and j=0. At 
these points fluxes are extrapolated according to the following third order 
formula. 

fj = 4/j4-i — 6 /j 4-2 + 4fj + 3 — /y-f 4 ( J — 0? 1) (3*5) 

fj + 1 = 4 fj — Qfj — i + 4/y — 2 — fj — 3 {j — J,J + 1) ( 3 * 6 ) 

When these extrapolations are applied to define the fluxes then the steps 
(3.3) and (3.4) are valid for all grid points j=l through J. 

This process is then applied to equation (2.8) with 

p — u 

} — U s p-{- p s U + up 

and 

h = -{U 8 p + p 8 u + U P)~£^ 

and to (2.9) with 

U — U 

f = U S U + — + c l{^Pf P s + ~^~2 (Pf 
and h — 0. 
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Artificial viscosity 


Similar to fluid dynamic calculations, in order to capture shock waves 
without oscillations, viscous damping terms are added to the difference 
equations. The previous numerical scheme described here is a dissipative 
scheme and therefore, we will show results with and without artificial 
viscosity. Von Neuman Richtmyer type viscous terms are used herein. If 
the equations (2.8) and (2.9) are written in the form: 

UL t + /(w)z =H{w) (3 .7) 

where w — 

then the artificial viscous term added to (3.7) is 



where V — 0(1). We differenced this quantity according to the following 
formula. 



This is a second order formula. This formula is used in both stages 
of our scheme up to the boundary. The added viscosity terms tend to 
reduce the accuracy of the scheme. Nevertheless, it provided better results 
than other types viscosity models we tried and gave a sharper shock with 
reasonably accurate solutions in the smooth regions. We shall see this later- 
in our comparison with the asymptotic method of reference 7. 
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4. Boundary Conditions. 


In this section we consider some questions concerning boundary con- 
ditions. First, we need a nonreflective boundary condition at the exit 
section of the duct. Next we need boundary conditions appropriate to 
the numerical scheme. These are accomplished by obtaining characteristic 
variables for the system (2.4) and (2.5). We linearize this system to get 

w t + Aw^ = 0 ( 4 . 1 ) 


where 



The eigenvalues of this matrix are 


X-j. = u s -{- c s and X_ = u s — c s 


We note that u s < c s thus X is strictly negative. The signs of these 
eigenvalues gives the characteristic directions of propagation flow. At x 
= 0, X_j_ > 0 gives the inflow direction and X_ < 0 gives the outflow 
direction. Similarly at x = L, X.}. > 0 gives the outflow and A— inflow 
directions respectively. The matrix formed by the eigenvectors is 


T = 



so that T X AT is diagonal. The characteristic variables are then 


( 4 . 2 ) 


then 


P , « 

Vi = — + - 

Ps C 8 

_ p u 

Ps Cs 


(4.4) 


Here Vi,V 2 corresponds to the eigenvalues X_|_ and X_ , respectively. At 
x = L V 2 is the inflow variable. Setting V 2 0, i.e. 

p U 

— =0 (4.5) 

Ps Cs 

is exactly the nonreflective condition stated by the general form in (2.13). 
In linear acoustics this is known as the impedence condition. 

For the numerical scheme it was found effective prescribing the boun- 
dary conditions in terms of these variables V\ and V 2 . At x = 0, l >2 is 

computed through an iteration. Let us call this value to be vZ, . Thus 


P_ _ U_ _ y c 
Ps C s 2 

But u is prescribed at x=0. 

U = f 

Solving (4.7) and (4.6) for u and p we have 

“ = /; p = p s {v 2 + 7 ) 


(4.6) 


(4.7) 


(4.8) 


Similarly, at x = L, we compute V\ . This gives 

p . u c 

— + — = 

Ps C G 


(4.10) 


We solve (4.5) and (4.10) to obtain the values of U and p at x = L. They 
are 

p = p s v\/2 
« — c e v\/2 
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Together with these conditions the solutions were started at a state 
of rest and iterated over 6 periods to obtain the results discussed in the 
next section. 

5. Discussion of results. 

The procedures developed in the previous sections were applied to a 
particular geometry called a Crocco-Tsien duct. A detailed description of 
the contour of the duct is available in reference 6. This contour is designed 
in such a way that, the mean flow accelerates linearly to Match number 
unity at the throat. In particular, for the examples given here the entry 
Mach number at the exit section was -.50 and at the throat -.90 
Here the “minus” sign denotes the flow in the negative x direction. 

Figure 1(a) shows a typical configuration of the duct. Figure 1(b) 
shows the area variation. This geometry has exit/throat ratio about 1.32, 
so that this area variation provides a gradual choking of the flow. In 
this case the Mach number distribution becomes as depicted in figure 
2. With this area variation and Mach number distribution, the steady 
flow equations satisfied by p s and u s (equations 2.5 and 2.6) were solved 
explicitly (see also [4]). 

As we discussed previously, the finite difference algorithm is compared 
with the asymptotic theory developed in reference 7. Since the typical 
nonlinear situation arises at higher sound pressure levels and Mach num- 
bers approaching unity, in this theory a small peaturbation parameter was 
chosen as (1— | Mt |), where Mt is the throat Mach number. Comparisons 
for a value for Mt = -.90 are given in figures 3 and 5 respectively. The 
strengths for an equivalent sound pressure source located at x=0 (figure 
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la) are roughly 149 and 156dB, respectively. Corresponding to equation 
(2.12) they have the form 

f(t) = A cos t, 

where A is the amplitude calculated according to the source strength. In 
both (149dB and 156dB) cases shocks were predicted in reference 7. In 
these figures the time history of the velocity over a period (2 7 r) is given. 
The velocity is normalized by (1 — | Mf |) 2 since the acoustic velocity 
is small in magnitude. Thus the actual value of the acoustic velocity 
is of the order 10 ® or less. This is exactly the reason a. higher order 
accurate scheme was necessary. The solid lines in these figures are the 
numerical solution and the other (see figures 3 and 5) are the asymptotic 
solutions computed at x = 0.75L and at the exit x = I> respectively. The 
finite difference calculations agree well with the asymptotic theory. The 
difference scheme we used itself is a dissipative scheme. We carried out 
the computations without artificial viscosity terms in the algorithm. For 
the case of 149dB source the results are shown in figure 7. The comparison 
is still good in the smooth regions. Results shown in figures 3 and 7 also 
validate the fact that the amount of added artificial viscosity did not affect 
the physics of the wave nature. 

Finally, figures 4 and 6 show the overall sound pressure level (dB) in 
the duct. Figure 4 corresponds to the shock case of 149dB sound source. 
This is a very weak shock case. At the exit we see a 2dB sound pressure 
level reduction. Figure 6 shows the sound pressure level for a 156dB source. 
In this case we see a sound pressure level reduction about 5dB at the exit. 
These results show that the higher the sound source level one obtains 
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substantial sound reduction through loss of energy due to shocks. 

Extension of this work in two dimensions is in progress by the authors. 
The results will be reported elsewhere. 
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